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We study the Bose-condensed ground states of bosons in a two-dimensional optical lattice in the presence 
of frustration due to an effective vector potential, for example, due to lattice rotation. We use a mapping to a 
large-5 frustrated magnet to study quantum fluctuations in the condensed state. Quantum effects are introduced 
by considering a 1/5 expansion around the classical ground state. The large-5 regime should be relevant 
to systems with many particles per site. As the system approaches the Mott insulating state, the hole density 
becomes small. Our large-5 results show that, even when the system is very dilute, the holes remain a (partially) 
condensed system. Moreover, the superfluid density is comparable to the condensate density. In other words, 
the large-5 regime does not display an instability to noncondensed phases. However, for cases with fewer 
than 1/3 flux quantum per lattice plaquette, we find that the fractional condensate depletion increases as the 
system approaches the Mott phase, giving rise to the possibility of a noncondensed state before the Mott phase 
is reached for systems with smaller 5 . 

PACS numbers: 03.75.Lm, 03.75.Mn, 75.10.Jm,75.10.-b, 75.45. +j 
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I. INTRODUCTION 

Bosonic atoms in optical lattices can display superfluid and 
Mott insulating phases. If the system is rotated, then, in the 
corotating frame, this is equivalent to introducing an effective 
magnetic field proportional to the rotation frequency'-^. This 
is not the only means to introduce a vector potential to a sys- 
tem of neutral atoms. This can also be achieved^^ through 
the interaction of atomic electric and magnetic moments with 
an external electromagnetic field (Aharonov-Casher and dif- 
ferential Aharonov-Bohm effects). For atoms trapped in an 
optical lattice in two distinct internal states, a scheme using 
two additional Raman lasers combined with the lattice accel- 
eration or inhomogeneous static electric field has also been 
proposed. 

Bosonic atoms in an optical lattice can be modeled by 
a Bose-Hubbard model. A vector potential introduces an 
Aharonov-Bohm phase for the boson hopping from site to site. 
The wave function is "frustrated" if the phase twists around 
each plaquette add up to 2na for some non-integer a. For a 
Bose condensate at a low effective magnetic field, this intro- 
duces vortices into the condensate. The presence of the optical 
lattice^'^ interferes with the formation of an Abrikosov vortex 
lattice and quantum fluctuations may be enhanced. Further, 
if the number of vortices becomes comparable to the num- 
ber of bosons, the system may enter into a fractional quantum 
Hall stateii2i2ii£iil. However, this requires a very high rotation 
frequency or a low atomic density which is hard to achieve 
experimentally. 

In this work, we will focus on the experimentally accessi- 
ble regime where a condensate still exists to examine whether 
there are any precursors to such states in a frustrated Bose 
condensate. We study a two-dimensional (2D) Bose-Hubbard 
model on a square lattice for a range of incommensurate fill- 
ing. In the regime of strong on-site interaction, the model is 
analogous to a quantum easy-plane ferromagnet and the frus- 
tration encourages spin twists, i.e., the formation of vortices 
in the ground state. We find the classical ground states using 
Monte Carlo methods and then we study the quantum fluctu- 



ations around the classical state. In other words, we work un- 
der the assumption that quantum effects do not change qual- 
itatively the nature of the ordering obtained for the classical 
ground states. Mathematically, this means that we will work 
in a large-5 generalization of the spin model and perform an 
expansion in IjS to obtain the quantum effects. Although our 
original model corresponds to small S , the large-S' approach 
canbejustified if the perturbative series in 1/5 converges 
In those cases, a spin wave calculation may give accurate re- 
sults. 

We will study how quantum fluctuations affect the order 
parameter, off-diagonal long-range order (ODLRO) and the 
superfluid fraction for different degrees of frustration for the 
whole range of incommensurate filling. In the spin analog, 
the incommensurate filling corresponds to a range of Zeeman 
field h up to some frustration-dependent critical field hda). 
Our calculations were made for a = 0, 1/4, 1/3, and 1/2. 

Our results show that the degree of Bose condensation de- 
creases as h increases toward he- However, it does not vanish 
at the limit of h - hda). This applies to several quantities 
that we have calculated: the reduction in the order parame- 
ter, the reduction in the largest eigenvalue of the density ma- 
trix, and the sum of the non-macroscopic eigenvalues of the 
density matrix. We also find similar conclusions for the su- 
perfluid fraction — frustration reduces the superfluid fraction 
in the comparison with the unfrustrated case but there is no 
vanishing of the superfluid fraction at any h < he- 

The paper is organized as follows. We will outline the 
model and the mapping to the quantum spin model in Sec. 
nil We describe the classical ground states (S — > oo) of the 
spin analog in Sec. [Ill] We introduce the excitations above 
the ground state in a l/S expansion in Sec. |IV] In Sees. IVl 
and lVIl we calculate the degree of condensation and superflu- 
idity in the system. We make conclusions about our study in 
the final section. 
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II. MODEL HAMILTONIAN 

For atoms trapped in a two-dimensional optical lattice, we 
can focus on a single-band lattice model if the tunneling f be- 
tween wells within the lattice is weak compared to the level 
spacings in each well. If the tunneling is also weak com- 
pared to the repulsive energy U for two atoms in one well, 
then strongly correlated ground states, such as the Mott insu- 
lator, appear as well as a superfluid state. 

Many different methods have been proposed to introduce 
frustration in the atomic motion. This can be done through 
rotating the system' or through the interaction of the atoms 
with an external electromagnetic field^"^. If there is only one 
species of bosonic atoms, then the system is described by a 
Bose-Hubbard model on a square lattice with a complex hop- 
ping matrix element: //nubbard - H^^^ + V with 

T = -f^(e*'fl]:a,+H.c.), (1) 

<'7> 

where ^ is the chemical potential and {ij) denotes nearest- 
neighbor sites / and j. The complex tunneling couplings ap- 
pear in the Hubbard Hamiltonian due to the presence of the ef- 
fective vector potential A. When an atom moves from a lattice 
site at R, to a neighboring site at Rj, it will gain an Aharonov- 
Bohm phase 




For neutral atoms with electric moments dg and a magnetic 
moments d„, in an external electromagnetic field {E, B), A - 
(dm + cteX S)lhc^^. For a rotating lattice, A - niQ. x r/h, 
where Q is the rotation frequency and m is the mass of the 
atom. In this work, we study the case of the uniform effective 
magnetic field B - V x A = Bz- Results will depend on the 
frustration parameter a, defined as the flux per plaquette in 
units of 2ji, 

plaq 

where the integration is over the surface of a lattice plaquette 
and the sum is performed anticlockwise over the edges of the 
square plaquette. This parameter is only meaningful between 
and 1 because a flux of In through a plaquette has no effect 
on the system. Frustration is maximal at a = 1 /2. 

In this paper, we will use a magnetic analogy as the frame- 
work to study the Bose-Hubbard problem. This is most easily 
motivated in the limit of Ujt — » 00, even though we will not be 
working directly in this limit. In such a limit, the site occupa- 
tion can be restricted to zero and one boson. Then, the Hilbert 
space of possible states can be mapped onto a spin-half XY 
model. The two S, states of the pseudospin correspond to 
whether a lattice contains a boson or not. 



The spin raising and lowering operators correspond to the 
creation and annihilation of hard-core bosons, respectively. 
This mapping is possible because hard-core bosons have the 
same commutation relations as 5 = 1/2 operators: operators 
on diff'erent sites commute but operators on the same site an- 
ticommute. The motion of the atoms translates to pseudospin 
exchange. The effective Hamiltonian is 

fi'^ff - -7 Z if^'-'StSj + H.c.) -hY,S) (4) 

<i;> J 

where J - 2t, Sf - S f + are the spin- 1/2 operators, and 
h - fi represents an effective Zeeman field. Note that this is a 
ferromagnet in the absence of frustration ((pij = 0). 

It is not simple to attack the infinite- 1/ limit of the problem 
of hard-core boson directly. Instead, we will relax the hard- 
core condition and allow for more than one boson on each site. 
We will allow 2S atoms on each site so that each site has 2S +1 
possible states. This corresponds to a spin-5 model with the 
Hamiltonian given in Eq. (|4]l. The relationship between the 
original bosons, a, and this spin-5 model is established via 
the Holstein-Primakoff representation: 

SI ^c'l(2S -c'lciy'\ 5t = c^Q-5. (5) 

where c,- are operators with bosonic commutations and are es- 
sentially the original bosons a, of the Bose-Hubbard model. 
The limit of 5 ^ 00 corresponds to the classical limit of the 
model. More specifically, we need 5 — > 00 while JS and h re- 
main constant so that exchange and Zeeman energies remain 
comparable. 

Mathematically, the large-5 limit provides a systematic 
way to control the quantum fluctuations in this problem. 
Quantum fluctuations can be introduced (see later) in a 1/S 
expansion under the assumption that those effects do not alter 
significantly the nature of the ordering obtained for the classi- 
cal ground states. We will present results to leading order in 
l/S (i.e., we do not set 5 = 1/2 afterward). Physically, the 
leading-order results in S should be relevant to optical lattices 
with many atoms per site on average. 

The relaxation of the maximum site occupancy to 25 from 
a model of hard-core bosons is not the only way to control 
correlations in the Bose-Hubbard model at weak tunneling. A 
similar methodology is to consider a dense but weakly inter- 
acting limit of the Bose-Hubbard model. With « being the av- 
erage boson density per site, this limit is given by f/ ^ and 
n ^ 00 while Uh remains constant"^. Then, one can develop a 
theory as an expansion in 1/h. This approach produces results 
very close to the 1/S expansion considered here. 

Note that our Hamiltonian has local gauge invariance. If we 
change the gauge, A — > A + Vx, then the Hamiltonian stays 
unchanged if the boson and spin operators pick up a phase 
change. 

<f>ij ^ e"^'-^'^<f>ij , hi e*fl,- , S; e'^'S; . (6) 

In the spin language, this corresponds to a rotation of;^^, in the 
xy plane in spin space. 
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Before proceeding to discuss the properties of this system, 
we point that we may generaHze this to an optical lattice con- 
taining two species of bosonic atoms, such as two hyperfine 
states. Let us denote the two species by cr -'[, j. This allows 
for more degrees of freedom in the model Hamiltonian. Two 
atomic species may, in general, see different lattice potentials 
so that the tunneling matrix elements and chemical potentials 
could be different for the two species. The Hubbard model for 
the two species would be of the form //Hubbard - //''"^ + T with 



i,(T 



(7) 



<r(ii) 



where the on-site interaction f/o-o-s the exchange interaction 
to-, the tunneling phase and the chemical potential jUo- have 
all acquired a dependence on the internal states of the bosons. 
If we specialize to the case of one atom per site with strong 
on-site interactions, we can rule out zero or double occupa- 
tion of each lattice site. In other words, the system should 
be a Mott insulator but the atom occupying each site can be 
of either internal state. Thus, each site has a spin-half de- 
gree of freedom: - a^^aii would create a t state and 



Sj - a[j^fl;T would create a X state. In this phase, the rela- 
tive motion of the two species of atoms is still possible: the 
motion of one species in one direction must be accompanied 
by the motion of the other species in the opposite direction. 
This counterflow keeps the occupation at one atom at each 
site. In the pseudospin language, this is simply spin exchange. 
Therefore, in this Mott phase for the overall density, we have 
again an easy-plane magnet. If we tune the interactions so 
that Uii = Uii = 2f/|^, then a perturbation theory in t/U 
brings us to the effective pseudospin Hamiltonian"*'^ described 
by Eq. © with J = A-t^tJU, h ^ 2 (ji^ - ^ii) + 8(f2 - f2)/f/, 

and cf>ij = (/>lj - cf>l. 

We can translate the phases of the single-species Hubbard 
model to this two-species system at unit filling. Superfluid- 
ity in the single-species Hamiltonian at an incommensurate 
filling corresponds to superfluidity for counterflow in the two- 
species problem at the commensurate filling of one atom per 
site but with different relative densities of the two species. 
The advantage of considering this two-species Mott insula- 
tor is that there may be more degrees of freedom in tuning the 
parameters of pseudospin Hamiltonian, including the explicit 
breaking of Sz — > -Sz spin symmetry. 



III. CLASSICAL GROUND STATES 

To determine the ground states of the pseudospin Hamilto- 
nian (|4|i, we consider first the 5 — > oo classical ground states 
for the spin system. We assume that h > without loss of 
generality. In the absence of the vector potential, the system 
is an easy-plane ferromagnet. For h < he - 4-JS , the ground 
state has a uniform magnetization in the xy plane in spin 
space. The xy component of the magnetization at each site 



is = [1 - (h/hc)^y^^. This xy magnetization corresponds 
to superfluidity in the original single-species Hubbard model. 
The z magnetization in the 5-' direction M; = N{S^} = Nh/hc 
corresponds to the number of atoms in the optical lattice mea- 
sured from half filling. For higher Zeeman fields [h > he), 
Mz becomes saturated and there is no xy magnetization: the 
lattice is a Mott insulator at one atom per site (or empty for 
h < -ha). 

In the presence of the vector potential, the ordering pattern 
of the classical ground state depends on the effective magnetic 
flux through each plaquette. This introduces vortices into the 
spin pattern. It also reduces the critical field h^ below which 
the xy magnetization is nonzero. As shown by Pazmandi and 
Domanski'^, he is given by is the maximal eigenvalue of the 
matrix JSe"^''. This is shown in Fig.[T] Note that this result 
for he is not restricted to the classical limit but applies for 
all values of the spin S . The spectrum of all the eigenvalues 
of this matrix as a function of the frustration parameter a is 
the Hofstadter spectrum.^" as discussed originally in terms of 
two-dimensional tight-binding electrons in the quantum Hall 
regime. 
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FIG. 1 : Critical value of the effective Zeeman field, hc(a), as a func- 
tion of the parameter a being the flux per plaquette in units of 2n. 
For h > hria) the lattice is a Mott insulator at one atom per site. 

Let us now turn to the classical ground states for h < he. 
Writing the local magnetization in spherical polars, (5^,) = 
S (sin 6i cos 0;, sin 6, sin 0;, cos 0,), the classical energy is given 
by: 

(ij) '■ 

(8) 

Minimizing this energy, we find that the ground-state values 
for (pi and 0,, <1), and 0,, must satisfy, for each site /, 

JS sin 0, ^ sin 0^ sin (o, - -i- 0^) = 

i=i+s 

JS cos 0, ^ sin 0/ cos (o, - <bj + 0^) - h sin 0, (9) 

j=i+S 

where the summation is taken over the four neighboring sites 
of /: j - i + S. The first equation conserves the spin cur- 
rent (or atomic current in the original Hubbard model) at each 



4 



node. The second specifies that there is no net effective Zee- 
man field causing precession around the z axis in spin space. 
In the original boson language, this ensures a uniform local 
chemical potential throughout the system (in the Hartree ap- 
proximation). The system has a local gauge invariance and we 
need to fix a gauge to perform our numerical calculations. We 
choose the Landau gauge A - B (0, x, 0) so that the Aharonov- 
Bohm phase (pij is zero on all horizontal bonds of the lattice. 
The classical ground states are obtained by using the Metropo- 
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FIG. 2: Ground state energy of the classical spin system as a function 
of the frustration parameter a (flux per plaquette divided by 2n) for 
difi'erent Zeeman fields h/JS = 0, 0.5, 1, 1.5, 2 and 2.5 (from top to 
bottom). The energy is symmetric around the point a = 1/2. 
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FIG. 3: Vortex patterns for (a) a = 1/3 and (b) o- = 1/2 (chequer- 
board configuration), with a being the flux per plaquette in units of 
2;r. For o' = 1/3 there are 2q = 6 degenerate states (vortices can be 
on three different 3x3 sublattices and along both diagonals). For 
a = 1/2 there are two degenerate states with vortices at one or the 
other diagonal. 



lis algorithm. For rational values of the frustration parameter 
or = p/q, the Monte Carlo simulations are done on nq x nq 
lattices with periodic boundary conditions. In most cases, we 
find that the periodicity of the ground state is q x q. How- 
ever, we also find ground states with the periodicity 2q x 2q 
in some cases. The ground-state energies as functions of the 
flux through a plaquette are shown in Fig.|2l 

We can also examine the vortex pattern in these ground 
states. The current on the bond joining sites / and j is given 
by: lij = (JS^/h)sm@iSm@jSm(<i>i - + (pij^. The circu- 
lation of these currents around each plaquette gives the vortex 
patterns. These are shown for a = 1/2,1 /3, and 1 /4 in Figs. 
Elandgl 

In case of a zero Zeeman field h - 0, the classical Hamilto- 
nian (|8]l has been studied extensively in the context of Joseph- 
son junction arrays in the presence of a perpendicular mag- 
netic field^'"^'^. Halsey^' showed that, for simple fractions in 
the range 1/3 < a < 1/2 (e.g., a = 1/2,1/3,2/5,3/7,3/8), 
the ground states have a constant current along diagonal stair- 
cases. Our results for h - Q agree with these previous stud- 
ies. For a general nonzero Zeeman field, the ground states we 
found for a - 112 and 1/3 also have currents in diagonal stair- 
cases. We cannot obtain analytic generalization of the Halsey 
solution for the case of finite h. We find the ground states by 
using the Metropolis algorithm. At finite h, the phase patterns 
for Of = 1 /2 and a - 1 /3 are similar to the phase patterns for 
the Halsey states aih - Q but has spatial variation around 
a finite average. 



FIG. 4: Vortex patterns for two ground states at a = 1/4 and h = Q. 
(a) Current pattern periodic on 4 x 4 square, phase pattern periodic 
on 8 X 8 square, (b) Current and phase patterns periodic on 4 x 4 
squares. 



The Halsey analysis does not cover cases when a < 1/3. At 
ff = 1 /4 and /i = we find two distinct ground state config- 
urations (Fig.lHl with the same energy in the agreement with 
previous results^Sr— . For both configurations, the current pat- 
terns are periodic on 4 x 4 square. However, the phase patterns 
do not have the same periodicity: it is 8 x8 periodic in the con- 
figuration shown in Fig.|4](a) but 4 x 4 in Fig.|4](b). We find 
states of the form [Fig.|4](b)] for general h when simulations 
are done on 4 x 4 lattices with periodic boundary conditions. 
Simulations done on larger An x An lattices at nonzero h give 
states that contain elements of both structures separated by 
domain walls. Similar results were found by Kasamatsu^i. 

IV. EXCITATION SPECTRUM 

In this section, we compute the excitations of the system 
using the spin-wave theory. Quantum effects are incorporated 
in the problem by considering finite values of S . We will per- 
form an expansion in powers of the parameter IjS and keep 
only the terms of the lowest order in 1 jS in the Hamiltonian. 
Even though we are interested in 5 ~ 0(1), the large-S' ap- 
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proach is in some cases justified due to the good convergence 
of the perturbative series Spin-wave approximation re- 
lies on an assumption that the introduction of the quantum 
fluctuations does not quahtatively change the nature of the or- 
dering obtained for classical ground state. We use this ap- 
proach to investigate whether the Bose condensate becomes 
unstable in any parameter regime. 

Starting from the classical ordered state, we use the 
Holstein-Primakoff transformation to represent the spin flips 
away from the classical ground state in terms of the bosonic 
operators. We will keep only the quadratic terms in the final 
bosonic Hamiltonian. It is convenient to introduce the oper- 
ators S, such that §>-J direction is parallel to the classical spin 
direction at each site 







sin 0, cos <!), 


sin 0; sin (5, 
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- sin O; 


cos O, 









- cos 0, cos 


- cos 0; sin O, 


sin 0, 



(10) 

and use the Holstein-Primakofi" representation of these new 
spin operators in terms of the bosonic operators, 

St = S!' + /Sr = (2S - P.biyi^bi, = 5 - b\bi. (11) 

Note that a gauge transformation corresponds to a rotation of 
the spin 5* around the z axis. Since these new spin variables 
are aligned with the classical spin configuration (whatever the 
choice of gauge), the new spin S is invariant under such rota- 
tion. Therefore, the bosonic operators, hi, are gauge invariant. 

Under assumption that the zero-point fluctuations are small 
so that the average number of spin flips at each site is small 
compared to 5 , we can approximate [ 1 -b'. hi / (25 )] ' as unity. 
The resulting Hamiltonian, to order 0(5'*'), is 

H - Ef'' + 2 [A.pihj - Athih] + H.c.) + 2 CiP.hi, (12) 

('7> i 

with 

Afj = ^ [(cos 0; COS &j + 1) Cij + /(cos 0/ + COS @j) Sjyj , 

Ci = 75 sin 0; ^ sin 0jC,y -H /i COS 0, (13) 

j=i+S 

where c,7 = cos(<l),-<l)j-i-0,y), Sij - sin(<l),-<l)j+0,y) and E^""" 
is the ground-state value of the classical energy [Eq. ©J. 
Note that all the coefficients in this Hamiltonian are gauge 
invariant, confirming our above conclusion that the bosonic 
operators, hi, are gauge invariant. 

This Hamiltonian also reduces correctly to the case of h > 
he (i.e., 0, - 0) when there is no need for realigning the axis of 
quantization [Eq. (fTOlll. In that case, the "anomalous" terms 
hh and Z?'^' in the Hamiltonian vanish. Then, the spin exci- 
tations are described by a tight-binding model with magnetic 
flux through the plaque ttes: 

Hh>H^, ^ -hNS -JSYj {e'''"'bih] + H.c.) -h /z ^ h]hi. (14) 

(ij) i 



This is diagonalized by the Hofstadter solution . The excita- 
tion spectrum has an energy gap of h - and the ground state 
corresponds to a vacuum of these excitations, i.e., there are no 
zero-point fluctuations in the ground state. 

For lower Zeeman fields {h < he), Hamiltonian (fT2] i con- 
taining the anomalous terms will have zero-point fluctuations 
which reduce the magnetization from the classical value. In 
the language of the original bosons, the fluctuations would 
deplete the condensate. The Hamiltonian can be diagonalized 
by a generalized Bogoliubov transformation, 

hi = {'^i'nO^m + Vl„al) , b] = {vimam + 

m m 

(15) 

for m - 1, . . ., / for a lattice of / sites. To ensure that the new 
operators a„, obey bosonic commutation relations, we require 
the matrices u and v to obey: uu' - vv' = 1 and uv^ - vu^ = 
0. To obtain a diagonalized Hamiltonian in terms of these 
new operators, we can write the part of the Hamiltonian ( fT2b 
quadratic in the bosonic operators as H = Mc, where M is 
a 2/ X 2/ matrix and c - (b,b') with b = (^1,^2, -■)■ Then, 
it can be shown that Hamiltonian ( fT2b is diagonalized into the 
form 

= £■() + ^ emfflam (16) 

m 

with eigenenergies 6,„ if we solve the eigenvalue problem, 

(m-^I,)^ = 0. (17) 

where q,„ = (ui,„, . . . ,un„„v\^^, . . . ,vl^J contains the co- 
efficients of the Bogoliubov transformation and Y., - 
{{1,0}, {0,-11). 

We computed the spectrum for 60 x 60, 120 x 120 and 
240 X 240 lattices with periodic boundary conditions, using 
the classical ground states from our Monte Carlo simulations 
discussed in the previous section. Our results for 60 x 60 lat- 
tices and the frustration parameters a - Q, 1/2, 1/3, and 1/4 
are shown in Fig.|5] Our result for a - 1/4 is calculated using 
the 4 x4 periodic classical ground state presented in Fig.|4jb). 

As can be seen in Fig. |5] at /1 < hda), the spectrum is 
gapless. The low-energy excitations are the Goldstone modes 
related to the spontaneous symmetry breaking of the global ro- 
tation symmetry in the xy-plane in spin space. In other words, 
the spin system has long-range magnetization in the xy plane 
in spin space. We can use (5^^) as the order parameter. In the 
language of the original bosonic model, this corresponds the 
breaking of U(l) symmetry due to Bose condensation. Above 
he, there is no symmetry breaking and we see an energy gap 
in the system proportional ioh- he as discussed above. 

The ground-state energy £0 [Eq. (fT6T l1 can be written as 
^■ciass _|. where IS.Eq = A -i- 2,„ e,„/2 is a quantum correc- 
tion to the classical ground-state energy [Eq. dHJ] with A = 
-75 cos((I)/ - + (pij) for /! = and -h 2; 1 /(2 cos 0,) 
for h Q. This quantum correction is of order 5 ^ while the 
classical energy is of order 5 and so the fractional change is 
small in the large-5 limit. We calculate the relative correc- 
tions A^o/^o""' for several lattice sizes (60 x 60, 120 x 120, 
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FIG. 5: Low energy excitation spectrum as a function of tlie Zeeman 
field h for 60 x 60 lattices with periodic boundary conditions for 
frustration o- = 0, 1/4, 1/3 and 1/2. Critical values are: hda = 
0) = 4, hAa = 1/4) = 2.828, h,(a = 1/3) = 2.732 and h,(a = 
1/2) = 2.828. Above h^, the spectrum has a finite energy gap. The 
spectrum is gapless for h < he indicating long-range order in the 
system. 



240 X 240) and extrapolate results to the thermodynamic limit 
shown in Fig. |6] As can be seen, the quantum correction de- 
creases to zero as the Zeeman field h approaches the critical 
value he- Above the ground state is the classical ground 
state containing no zero-point fluctuations. 
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FIG. 6: Quantum correction to the ground-state energy as a function 
of h/hc(a) for a = 0, 1/2, 1/3 and 1/4. hdcr) is the critical value of 
the Zeeman field h for a given frustration parameter a. 



our boson problem. Since we are considering a lattice system 
above half filling, it is more meaningful to consider the con- 
densation of vacancies because this is the most appropriate 
description as h approaches h^. (For the two-species model 
with counterflow superfluidity, we are considering the con- 
densation of the minority species.) The hole density matrix is 
defined as p^. - {aja\). The existence of a macroscopic eigen- 
value, /Vo, coiTesponds to Bose-Einstein condensation. The 
sum of aU non-macroscopic eigenvalues gives the number of 
holes not in condensate and we can define the fractional con- 
densate depletion as the ratio of the non-macroscopic sum to 
the total number of holes /Vh which is the trace of the density 
matrix. 

In the analog of the easy-plane magnet, we should study the 
spin-spin correlation function for the spin components in the 
xy plane: pj, = {SjS^j). ODLRO corresponds to a non-zero 
xy magnetization which is the analog of Bose condensation. 
In the large-5 limit, Pji/2S is the analog of the bosonic hole 
density matrix p^. for h close to he- 

The macroscopic eigenvalue for our spin-spin correlation 
function is, to the leading order in S , given by the classical 
value /Vg'"'*^ = E/('"^'' )^' where m^- is the classical value of the 
magnetization at site We present below our results for con- 
densate and the depletion of the condensate, i.e., zero-point 
fluctuations which decrease the magnetization in the ground 
state. 

The above discussion needs to be modified in the presence 
of a vector potential because the density matrices, p and p^, 



are not gauge-invariant quantities: pji e 



'■^'^Pji under the 



gauge transformation [Eq. (|6])]. However, we can construct 
gauge-invariant analogs. Moreover, the eigenvalues of p and 
p^ are gauge invariant even though the corresponding eigen- 
vectors are not. Consider first the spin-spin correlation func- 
tion in the ground state 

pji = {SiSp^pf- + 6pji 
pf^' ^ with i/r; = 5e'*' sin0; (18) 

where p^'^^'* is the classical value of the density matrix (of or- 
der S^) and t//i is the classical value of the order parameter 
(of order S) (Sf). The order parameter itself is reduced by 
quantum fluctuations, 

{§!) = ^Kl - A,) , A,- = 1 2 Iv^J. (19) 



The correction Sp to the density matrix is given by: 



Spji - -p;.|-^(A,. + A,.) + L-C^^-'^O ^ (20) 



V. DENSITY MATRIX 

In this section, we will examine ODLRO in the density 
matrix^Si^^. Consider first the case without a vector poten- 
tial. A macroscopically large eigenvalue of the density ma- 
trix Pji signals the existence of Bose-Einstein condensation for 



where = m,„ + v',„ + cos 0,(v,„ - m/„), with m,>, and v',„ being 
the coefficients for the Bogoliubov transformation [Eq. ( fTSl) ]. 
This density matrix is not invariant under a gauge transforma- 
tion. We obtain a gauge-invariant version of the density ma- 
trix by expressing it with respect to a gauge-covariant basis. 
The most natural basis is the basis formed by the eigenvectors 
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of the classical density matrix p'^'^^^ The eigenvector corre- 
sponding to the largest eigenvalue is simply i^,, 

Y,plr^i = A^s'^^Vi with 



where A^g'"^^ is simply the classical value of the sum of the 
square of the xy magnetization on each site. It is on the 
order of NS ^ at /i = and tends to zero as h reaches he- All the 
other eigevectors of p'^'"^^ have eigenvalues of zero. Using an 
orthonormal set of these eigenvectors as columns for a unitary 
matrix U, we can construct a unitary transformation for the 
density matrix {p p, etc.), 



p = U'ipU = p'^'"' + dp . 



(22) 



where p'^'^^^ - diag(A^()'^^% 0, . . . , 0). Under the gauge transfor- 
mation [Eq. (|6])], all the eigenvectors of pji pick up a phase 
change, e.g., t/zj — > e^'^'t/zj so that U/j — > e^'^'Uij. It is easy to 
check that this compensates for the phase change in py, so that 
pij — * Pij- Consequently, all the quantities obtained from the 
matrix p are gauge-invariant and therefore physically mean- 
ingful. In this section, we calculate the effect of quantum fluc- 
tuations on the density matrix. This requires only the eigen- 
values of p. They are in fact the same as the eigenvalues of 
p because the two density matrices are related by a unitary 
transformation. 

We will now present our numerical results. First of all, we 
present the classical solution for the number of atoms in the 
condensate, N^^'^^'^^, as given by Eq. ( 1211 1. This is shown in 
Fig. Ill We see that this decreases to zero as h is increased to 
hcia). 
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FIG. 7: The condensate density per site, Ng'"^^"^'^' / 1 , in the classical 
limit for a = 0, 1/4, 1/3, 1/2. (/ = number of lattice sites.) hda) is 
the critical value of the Zeeman field h for a given frustration param- 
eter a. 

Next, we compute the quantum corrections to the classi- 
cal solution. In the large-S' expansion, these corrections are 
small and the leading corrections are of order 1/5 compared 
to the classical limit. We have computed this leading-order 



correction and present results in terms of the correction to the 
classical limits as fractions of the classical solution. 

We can exploit the large-5 expansion to compute the eigen- 
values of the density matrix. We start with calculating the 
quantum correction to the non-degenerate macroscopic eigen- 
value, A^(). Since p^j"^^ is larger than Spji by an order in S , 
we can calculate the eigenvalues of p by treating Sp in pertur- 
bation theory. The first-order correction to A^o is then given 
by 

^^o^-^X-^^^^'V-A.-^PH (23) 

if the first basis vector for dp is chosen to be the one corre- 
sponding to the classical solution i/r. This correction is of order 
S , as opposed to order for the classical value. Our results 
for AA^o as a fraction of A^^'^^'* are shown in Fig. [8] We see that 
the reduction in A^o is largest at /; = and decreases to zero at 
the critical fields hda). The vanishing of quantum corrections 
as h he (0; — » 0) can be seen directly from the coefficients 
A" of the anomalous terms in Hamiltonian (fT2] ) which are re- 
sponsible for the zero-point fluctuations in the ground state. 
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FIG. 8: Quantum correction AA'o to the the macroscopic eigenvalue 
of the density matrix as a function of hi hda) for a = Q, 1/4, 1/3 
and 1 12. Results have been extrapolated to the thermodynamic limit 
(L — » oo). hc(a) is the critical value of the Zeeman field h for a given 
frustration parameter a. 

We can also calculate the sum of the non-macroscopic 
eigenvalues, A^out- This corresponds to the condensate deple- 
tion in the original boson problem. In the 5 — > oo limit for 
a lattice with / sites, the / - 1 non-macroscopic eigenvalues 
are ah zero. The first-order quantum corrections can be ob- 
tained using degenerate perturbation theory — we can obtain 
the eigenvalues as the eigenvalues of the (/ - l)-dimensional 
submatrix 6pji for i,] - 2, . . . , / which excludes the macro- 
scopically occupied state. The sum of these eigenvalues is 
simply the trace of the submatrix: 



A^o, 



(24) 



!#1 



Again, A^out °^ S is one order smaller in S than N^'^^'^. We find 
that, just as classical condensate density (N^^^^'^/I) vanishes as 
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h hcia), the out-of-condensate number, A^out, also vanishes 
as /z ^ hcia). However, the ratio of the two quantities remains 
finite. This ratio, A^out/A^(f^^% is fractional depletion of the 
condensate. This quantity is one of interest in experiments 
which measure the degree of Bose-Einstein condensation by 
observing the time of flight of expanding condensates. Our 
resuhs for this fractional depletion A^out/A^o'^^\ rescaled by 5, 
are shown in Fig.|9] The occupation of these non-macroscopic 
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FIG. 9: Fractional depletion K^JN^'''"' for a = 0, 1/4, 1/3 and 1/2 
and as a function of h/hda). hc(a) is the critical value of the Zee- 
man field h for a given frustration parameter a. Results have been 
extrapolated to the thermodynamic limit (L — » oo). 



modes is also due to the anomalous terms in the Hamiltonian. 
This again should vanish as h — > he- However, Fig. |9] shows 
that the occupation remains a. finite fraction of N^^^^ even at 
the critical field he- In terms of the original boson model, 
this result suggests that condensate depletion remains a finite 
fraction of the total number of holes even as the hole density 
decreases to zero at he- Our results at zero frustration agrees 
with previous work'^-^^. 

We observe that this fractional depletion decreases monot- 
ically as we increase the Zeeman field h from zero to he for 
a = and 1/2. For a - 1/3, the fractional depletion appears 
to have zero slope as a function of h near he- Interestingly, 
for a - 1/4, the relative depletion becomes a non-monotonic 
function of the Zeeman field — the fractional depletion in- 
creases when he is approached. In fact, if we formally set 
5 = 1/2, the condensate depletion even reaches unity before 
h reaches he- As we will see in the next section, this change 
in behavior for a - 1 /4 is also seen in the superfluid fraction. 
We discuss this further in our concluding remarks. 

We note that A^out -AN[). In other words, the trace of 
the density matrix changes due to quantum fluctuations. This 
means that, in the quantum magnet, there is more than one 
possible measure of "condensation" in the ground state. The 
discrepancy can be traced to the quantum fluctuations for S~ 
at each site: Trp = ZK^+^r) = ^.[^(^ + l)-<(5;)2> + <5;)]. 
For 5 = 1/2, this is simply 2,(1/2 + {S'j}), corresponding to 
the total boson number in the original model which is a con- 
served quantity. However, for any S > 1/2, the mean-square 
fluctuation in the local z-component will alter the total trace 
of the density matrix. In other words, this is an artifact of our 



large-5 generalization of the model. In the above, we have 
compared A^out with the macroscopic eigenvalue A^o - A^q'^^^ 
Strictly speaking, in order to discuss the depletion of the hole 
condensate in the original boson model, we should use the 
analogue for the hole density matrix and then divide the num- 
ber of holes in the system. As discussed above, the correspon- 
dence is simple near /i^: we should consider Nout/2S com- 
pared to 2,(5 - (5?» = S 2i(l - cos 0,). This is quahtatively 
similar to the results plotted in Fig.|9] 



VI. SUPERFLUID DENSITY 

Bose-Einstein condensation can be defined in equilibrium. 
On the other hand, superfluidity is related to the transport 
properties of the system. Those two phenomena are related 
through the phase of the macroscopic wave function (order 
parameter). The superflow occurs when the phase of the wave 
function varies in space. In this section, we calculate the su- 
perfluid density for our system as a response to an external 
phase twist. The superfluid density, a characteristic quantity 
that describes the superfluid, measures the phase stiffness un- 
der an imposed phase variation and differs from zero only in 
the presence of the phase ordering. We find the superfluid 
fraction following the calculations of Roth and Burnett^^ and 
Rey et al:^ where the superfluid density is calculated for the 
Bose-Hubbard model with real couplings. Our results show 
that the superfluid fraction is reduced in the presence of the 
frustration. 

The superfluid density introduced by considering a change 
in the free energy of the system under imposed phase 
variations'**"^" is equivalent to the helicity modulus-'" which 
diff'ers from zero only for ordered-phase configurations and 
is consequently an indicator of the long-range phase coher- 
ence of the system. The definition is also equivalent to the 
definition of the superfluid density in terms of the wind- 
ing numbers which is used in the path-integral Monte Carlo 
methods^ '"■'^ and to Drude weight or charge stiffness which 
describes d.c. conductivity ■'''"■'^. 

Let us consider a system of size Lj in the x direction. One 
way to achieve the phase twist is to impose the twisted bound- 
ary conditions on the wave function describing the system. If 
we assume that the phase twist is imposed along the x direc- 
tion the twisted boundary conditions are 

{A,..., fi + LA ...) = e'^^'^ {A,..., A, ...) (25) 

with respect to all coordinates of the wave function. Let us 
introduce a unitary transformation 



with ^^xir + L^ -X{r)- 



(26) 



The untwisted wave function which satisfies the periodic 
boundary conditions ^(ri, + Li.£, ...) - ^^(ri, r,, ...) 
is related to the twisted wave function via the unitary trans- 
formation as I*!**) = U^'^). The Schrodinger equation 
for the system with twisted boundary conditions, H\^'^) - 
E''\^'^), can then be rewritten as H^\^') = E^\^') where the 
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twisted Hamiltonian is 



(27) 



In other words, the eigenvalues of the twisted Hamiltonian 
with periodic boundary conditions are the same as eigenvalues 
of the original Hamiltonian with twisted boundary conditions. 
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FIG. 10: Superfluid density as a fraction of the classical condensate 
density Nfl'"^^/I as a function h/h^(a) for the frustration parameter 
a = 0, 1/4, 1/3 and 1/2. hc(a) is the critical value of the Zeeman 
field h for a given frustration parameter a. 

The superfluid velocity is proportional to the order- 
parameter phase gradient and an additional phase variation 
X{r) will change the superfluid velocity by Av^ = hVx{f)lm 
in the continuous system. When the imposed phase gradi- 
ent is small so that other excitations except increase in the 
velocity of the superflow can be neglected the change in the 
ground-state energy can be approximated by AEg - -P- Avs + 
MsiAvs)^/2, with Ms - inNs being the total mass of the super- 
fluid part of the system. Here we choose a linear phase vari- 
ation along the x direction, ;i'(r) - ix/L^. Replacing tp-jlm 
for the continuous system by 7/2 for our 2D discrete lattice we 
obtain the following expression for the superfluid density^ 



0=0 



(28) 



where /^_y - Lv,y/« with a being the lattice spacing. The 
twisted Hamiltonian is of the same form as the untwisted one 
only with replaced by (pjj - O. Under assumption that the 
phase twist <«c ;r we can calculate the ground state energy 
of the twisted Hamiltonian perturbatively. Expanding e'*^^' 
up to the second order in O the twisted spin Hamiltonian be- 



= 5/1-//1- - C)2f 1-/2/2. Calculating the ground-state en- 
ergy for the system with imposed small twist within the sec- 
ond order perturbation theory and using Eq. (l28T l, we obtain 
the following expression for the superfluid density as a frac- 
tion of the condensate density, /j - IJynJNo: 



1 



v#0 



(30) 

where A^o - Nf'^^^ in the large-5 Umit and i/Zy are eigenstates of 
original untwisted Hamiltonian with v = labeling the ground 
state. In terms of the original boson model, A^o corresponds 
to the number of condensed particles or holes (for h < Q or 
h > 0). The first term corresponds to the diamagnetic response 
of the condensate while the second term corresponds to the 
paramagnetic response involving excited states. 

The results obtained for the superfluid fraction within the 
Bogoliubov approximation are shown in Fig. [TOl The lead- 
ing term due to quantum effects comes from the paramagnetic 
term in Eq. dSOl ). This is of order 5 In the absence of frustra- 
tion (a = 0), the system is homogenous and the system con- 
serves momentum. This means that the eigenstates are Bloch 
states corresponding to different momenta. As a result, the 
current matrix element in Eq. ( |30] |, which cannot couple dif- 
ferent momenta, vanishes. Moreover, the kinetic energy in the 
ground state is in itself proportional to A^o- In the boson model, 
this means that the superfluid fraction corresponds simply to 
the kinetic energy per hole. This is a quantity which is in- 
dependent of h and so the superfluid density is the same as 
the condensate density in the laige-S limit at zero frustration. 
(However, 1 /S corrections will change the result, giving a su- 
perfluid density larger than the condensate density for general 
h, but /s — » 1 as /; — > he-) Similarly, the current matrix ele- 
ment vanishes for the fully frustrated case (a - 1/2). In this 
case, frustration reduces the superfluid fraction in a - 1/2 
case to around 70%. For a - 1/3 and 1/4, an increase in 
the Zeeman field h results in a larger reduction in the fraction 
/, at values of h closer to hda). That can be seen in Fig.fTOl 
for the inhomogeneous cases of a = 1/3 and 1 /4. As for the 
condensate depletion, we note that the superfluid density as a 
fraction of the condensate density does not vanish as h he- 

We also note that the superfluid density behaves differently 
for a = 1/3 and 1/4 compared to a = and 1/2. The same 
qualitative change in behavior was observed for the conden- 
sate depletion calculated in Sec. [V] 



VII. CONCLUSION 



/v 2/2 " 



(29) 



where - iJ Tjii^"^"*'SfS ^^^^ - H.c.)/2 is the paramagnetic 
current operator and T^- - -J + H.c.)/2 cor- 

responds to the kinetic-energy operator for the hopping in the 
X direction. The terms in the Hamiltonian above that con- 
tain the twist angle can be treated as a small perturbation 



We have studied the ground state for bosonic atoms in a 
frustrated optical lattice by mapping the problem to a frus- 
trated easy-plane magnet. Using a large-5 approach, we 
further introduce quantum effects under the assumption that 
those effects do not change qualitatively the nature of the or- 
dering obtained for the classical ground states. We examined 
our results for any precursor to the non-superfluid or uncon- 
densed states. 
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We have found that frustration can decrease the depletion of 
the condensate and the superfluid fraction. However, the frac- 
tional depletion of the condensate and the superfluid fraction 
remain finite for all incommensurate filhng [h < hda)]. The 
behavior of the fractional condensate depletion and superfluid 
fraction as a function of filling has interesting behavior. We 
find that the cases of a = and 1/2 behave differently from 
the cases of o- = 1/3 and 1/4. Surprisingly, for the cases of 
smaller a, the fractional condensate depletion becomes a non- 
monotonic function of the filhng, decreasing as we increase 
h from zero but eventually increases as h —> h^. In fact, if 
we formally set 5 = 1 /2, then the computed fractional deple- 
tion exceeds 100% for the a = 1/4 case as h approaches he- 
We also have some evidence that the same behavior occurs 
in the a = 1/6 case for small system sizes. In other words, 
our results raise the possibility, for a < 1 /4, of a second-order 
phase transition to a non-condensed state where quantum fluc- 
tuations are large enough to destroy Bose condensation. It is 



intriguing to note that this case does not have a Halsey-type 
classical ground state and in fact has two degenerate ground 
states with different phase patterns. One can speculate that the 
motion of domain walls between the two different phase pat- 
terns may contribute to a route to decondensation and/or loss 
of superfluidity. 

Finally, we note that fractional quantum Hall states are ex- 
pected when the number of vortices becomes comparable to 
the number of atoms or holes in the Bose-Hubbard model. In 
our large-5 theory, the boson number is proportional to S and 
so the quantum Hall regime, if it exists in such a theory, exists 
only when h - he ~ l/S . Therefore, one might expect the 
condensate depletion or the reduction in the superfluid frac- 
tion to be large as h ^ he. We do not find this directly in our 
perturbative theory in 1/S . However, our results for the fluc- 
tuations around non-Halsey-type ground states suggest that an 
instabihty to a non-condensed state may be possible. 
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